#
# analyzeburg.r
# Replication and reanalysis after multiple imputation of Burgoon dataset, JCR 2006.
#   This code creates all the quantities of interest to compare the effects of 
#   imputation, as well as the graph included in the paper as figure 6.
#
# 23/2/10, jH
#

library(MASS)
library(Zelig)

model.burgoon<-function(rd){

  rd$year<-as.factor(rd$year)

  # Compute Coefficients

  output1a<-glm.nb(terrorinclead~totspendrevlog + govleft + democ + poplog + govcap + conflict + tradelog 
                          + europe +africa + asia + america + terrorinc + year, data=rd)

  output2a<-glm.nb(terrorinclead~transferslog + govleft + democ + poplog + govcap + conflict + tradelog 
                          + europe +africa + asia + america + terrorinc + year, data=rd)

  output3a<-glm.nb(terrorinclead~welfarelog + govleft + democ + poplog + govcap + conflict + tradelog 
                          + europe +africa + asia + america + terrorinc + year, data=rd)

  output1b<-glm.nb(terrorinclead~totspendrevlog + govleft + democ + poplog + govcap + conflict + tradelog 
                          + europe +africa + asia + america, data=rd)

  output2b<-glm.nb(terrorinclead~transferslog + govleft + democ + poplog + govcap + conflict + tradelog 
                          + europe +africa + asia + america, data=rd)

  output3b<-glm.nb(terrorinclead~welfarelog + govleft + democ + poplog + govcap + conflict + tradelog 
                          + europe +africa + asia + america, data=rd)

  # Compute Standard Errors

  coef<-cbind(output1a$coefficients[1:13],output2a$coefficients[1:13],output3a$coefficients[1:13],c(output1b$coefficients[1:12],NA),c(output2b$coefficients[1:12],NA),c(output3b$coefficients[1:12],NA))
  coef<-rbind(coef[2:nrow(coef),],coef[1,])

  se1a<-summary(output1a)

  se2a<-summary(output2a)

  se3a<-summary(output3a)

  se1b<-summary(output1b)

  se2b<-summary(output2b)

  se3b<-summary(output3b)

  se<-cbind(se1a$coefficients[1:13,2],se2a$coefficients[1:13,2],se3a$coefficients[1:13,2],c(se1b$coefficients[1:12,2],NA),c(se2b$coefficients[1:12,2],NA),c(se3b$coefficients[1:12,2],NA))
  se<-rbind(se[2:nrow(se),],se[1,])

  return(list(coef,se))
}


# simple function hardcoded for m=10 imputations to combine results by Rubin's rules.
my.mi2<-function(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,s1,s2,s3,s4,s5,s6,s7,s8,s9,s10){
  meanb<-(b1+b2+b3+b4+b5+b6+b7+b8+b9+b10)/5
  s.d<-((b1-meanb)^2+ (b2-meanb)^2+ (b3-meanb)^2+ (b4-meanb)^2+ (b5-meanb)^2+ (b6-meanb)^2+ (b7-meanb)^2+ (b8-meanb)^2+ (b9-meanb)^2+ (b10-meanb)^2)/(9)
  s.m<-(s1^2+ s2^2+ s3^2+ s4^2+ s5^2+ s6^2+ s7^2+ s8^2+ s9^2+ s10^2)/10
  seb<-(s.m + s.d*(1.1))^0.5
  return(list(meanb,seb))
}


x<-load("burgoonsubset.RData")

lwd<-model.burgoon(reasonabledata)

imp1<-read.csv("p3ny1.csv")
imp2<-read.csv("p3ny2.csv")
imp3<-read.csv("p3ny3.csv")
imp4<-read.csv("p3ny4.csv")
imp5<-read.csv("p3ny5.csv")
imp6<-read.csv("p3ny6.csv")
imp7<-read.csv("p3ny7.csv")
imp8<-read.csv("p3ny8.csv")
imp9<-read.csv("p3ny9.csv")
imp10<-read.csv("p3ny10.csv")

p1<-model.burgoon(imp1)
p2<-model.burgoon(imp2)
p3<-model.burgoon(imp3)
p4<-model.burgoon(imp4)
p5<-model.burgoon(imp5)
p6<-model.burgoon(imp6)
p7<-model.burgoon(imp7)
p8<-model.burgoon(imp8)
p9<-model.burgoon(imp9)
p10<-model.burgoon(imp10)


imp<-my.mi2(p1[[1]],p2[[1]],p3[[1]],p4[[1]],p5[[1]],p6[[1]],p7[[1]],p8[[1]],p9[[1]],p10[[1]],p1[[2]],p2[[2]],p3[[2]],p4[[2]],p5[[2]],p6[[2]],p7[[2]],p8[[2]],p9[[2]],p10[[2]])

cat("Original Parameters From Replication\n")
print(lwd[[1]],digits=2)

cat("Parameters After Imputation\n")
print(imp[[1]],digits=2)

plot(imp[[1]][1:6,],lwd[[1]][1:6,],xlim=c(-1,1),ylim=c(-1,1),xlab="imputations",ylab="listwise del.")
par(new=TRUE)
plot(imp[[1]][7,],lwd[[1]][7,],xlim=c(-1,1),ylim=c(-1,1),xlab="imputations",ylab="listwise del.",col="red")
par(new=TRUE)

imp.ci.lower=imp[[1]]-1.96*imp[[2]]
imp.ci.upper=imp[[1]]+1.96*imp[[2]]

lwd.ci.lower=lwd[[1]]-1.96*lwd[[2]]
lwd.ci.upper=lwd[[1]]+1.96*lwd[[2]]

for(i in 1:6){
for(j in 1:6){
  par(new=TRUE)
  plot(rbind(imp[[1]][i,j],imp[[1]][i,j]),rbind(lwd.ci.lower[i,j],lwd.ci.upper[i,j]),xlim=c(-1,1),ylim=c(-1,1),xlab="imputations",ylab="listwise del.",type="l")
  par(new=TRUE)
  plot(rbind(imp.ci.lower[i,j],imp.ci.upper[i,j]),rbind(lwd[[1]][i,j],lwd[[1]][i,j]),xlim=c(-1,1),ylim=c(-1,1),xlab="imputations",ylab="listwise del.",type="l")
}
}


for(i in 7){
for(j in 1:6){
  par(new=TRUE)
  plot(rbind(imp[[1]][i,j],imp[[1]][i,j]),rbind(lwd.ci.lower[i,j],lwd.ci.upper[i,j]),xlim=c(-1,1),ylim=c(-1,1),xlab="imputations",ylab="listwise del.",type="l",col="red")
  par(new=TRUE)
  plot(rbind(imp.ci.lower[i,j],imp.ci.upper[i,j]),rbind(lwd[[1]][i,j],lwd[[1]][i,j]),xlim=c(-1,1),ylim=c(-1,1),xlab="imputations",ylab="listwise del.",type="l",col="red")
}
}


abline(a=0,b=1)
abline(h=0)
abline(v=0)


m<-(1:6)*0
mlwd<-(1:6)*0

# Model 1

z.out <- zelig(terrorinclead~totspendrevlog + govleft + democ + poplog + govcap + conflict + tradelog 
           + europe +africa + asia + america + terrorinc + year, model="negbin", 
           data=mi(imp1,imp2,imp3,imp4,imp5,imp6,imp7,imp8,imp9,imp10))

z.lwd.out <- zelig(terrorinclead~totspendrevlog + govleft + democ + poplog + govcap + conflict + tradelog 
           + europe +africa + asia + america + terrorinc + year, model="negbin", 
           data=reasonabledata)

x.low <- setx(z.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)-sd(reasonabledata$tradelog,na.rm=TRUE))
x.high <- setx(z.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)+sd(reasonabledata$tradelog,na.rm=TRUE))

x.lwd.low <- setx(z.lwd.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)-sd(reasonabledata$tradelog,na.rm=TRUE))
x.lwd.high <- setx(z.lwd.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)+sd(reasonabledata$tradelog,na.rm=TRUE))

s.out <- sim(z.out,x=x.low,x1=x.high)
s.lwd.out <- sim(z.lwd.out,x=x.lwd.low,x1=x.lwd.high)

junk<-sort(s.out$qi$fd)
ci1<-c(junk[round(length(junk)*0.05)],junk[round(length(junk)*0.95)])
m[1]<-mean(junk)

junk<-sort(s.lwd.out$qi$fd)
cilwd1<-c(junk[round(length(junk)*0.05)],junk[round(length(junk)*0.95)])
mlwd[1]<-mean(junk)

# Model 2

z.out <- zelig(terrorinclead~transferslog + govleft + democ + poplog + govcap + conflict + tradelog 
                          + europe +africa + asia + america + terrorinc + year, model="negbin", 
           data=mi(imp1,imp2,imp3,imp4,imp5,imp6,imp7,imp8,imp9,imp10))
z.lwd.out <- zelig(terrorinclead~transferslog + govleft + democ + poplog + govcap + conflict + tradelog 
                          + europe +africa + asia + america + terrorinc + year, model="negbin", 
           data=reasonabledata)

x.low <- setx(z.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)-sd(reasonabledata$tradelog,na.rm=TRUE))
x.high <- setx(z.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)+sd(reasonabledata$tradelog,na.rm=TRUE))

x.lwd.low <- setx(z.lwd.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)-sd(reasonabledata$tradelog,na.rm=TRUE))
x.lwd.high <- setx(z.lwd.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)+sd(reasonabledata$tradelog,na.rm=TRUE))

s.out <- sim(z.out,x=x.low,x1=x.high)
s.lwd.out <- sim(z.lwd.out,x=x.lwd.low,x1=x.lwd.high)

junk<-sort(s.out$qi$fd)
ci2<-c(junk[round(length(junk)*0.05)],junk[round(length(junk)*0.95)])
m[2]<-mean(junk)

junk<-sort(s.lwd.out$qi$fd)
cilwd2<-c(junk[round(length(junk)*0.05)],junk[round(length(junk)*0.95)])
mlwd[2]<-mean(junk)

# Model 3

z.out <- zelig(terrorinclead~welfarelog + govleft + democ + poplog + govcap + conflict + tradelog 
                          + europe +africa + asia + america + terrorinc + year, model="negbin", 
           data=mi(imp1,imp2,imp3,imp4,imp5,imp6,imp7,imp8,imp9,imp10))
z.lwd.out <- zelig(terrorinclead~welfarelog + govleft + democ + poplog + govcap + conflict + tradelog 
                          + europe +africa + asia + america + terrorinc + year, model="negbin", 
           data=reasonabledata)

x.low <- setx(z.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)-sd(reasonabledata$tradelog,na.rm=TRUE))
x.high <- setx(z.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)+sd(reasonabledata$tradelog,na.rm=TRUE))

x.lwd.low <- setx(z.lwd.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)-sd(reasonabledata$tradelog,na.rm=TRUE))
x.lwd.high <- setx(z.lwd.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)+sd(reasonabledata$tradelog,na.rm=TRUE))

s.out <- sim(z.out,x=x.low,x1=x.high)
s.lwd.out <- sim(z.lwd.out,x=x.lwd.low,x1=x.lwd.high)

junk<-sort(s.out$qi$fd)
ci3<-c(junk[round(length(junk)*0.05)],junk[round(length(junk)*0.95)])
m[3]<-mean(junk)

junk<-sort(s.lwd.out$qi$fd)
cilwd3<-c(junk[round(length(junk)*0.05)],junk[round(length(junk)*0.95)])
mlwd[3]<-mean(junk)

# Model 4

z.out <- zelig(terrorinclead~totspendrevlog + govleft + democ + poplog + govcap + conflict + tradelog 
                          + europe +africa + asia + america, model="negbin", 
           data=mi(imp1,imp2,imp3,imp4,imp5,imp6,imp7,imp8,imp9,imp10))
z.lwd.out <- zelig(terrorinclead~totspendrevlog + govleft + democ + poplog + govcap + conflict + tradelog 
                          + europe +africa + asia + america, model="negbin", 
           data=reasonabledata)

x.low <- setx(z.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)-sd(reasonabledata$tradelog,na.rm=TRUE))
x.high <- setx(z.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)+sd(reasonabledata$tradelog,na.rm=TRUE))

x.lwd.low <- setx(z.lwd.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)-sd(reasonabledata$tradelog,na.rm=TRUE))
x.lwd.high <- setx(z.lwd.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)+sd(reasonabledata$tradelog,na.rm=TRUE))

s.out <- sim(z.out,x=x.low,x1=x.high)
s.lwd.out <- sim(z.lwd.out,x=x.lwd.low,x1=x.lwd.high)

junk<-sort(s.out$qi$fd)
ci4<-c(junk[round(length(junk)*0.05)],junk[round(length(junk)*0.95)])
m[4]<-mean(junk)

junk<-sort(s.lwd.out$qi$fd)
cilwd4<-c(junk[round(length(junk)*0.05)],junk[round(length(junk)*0.95)])
mlwd[4]<-mean(junk)

# Model 5

z.out <- zelig(terrorinclead~transferslog + govleft + democ + poplog + govcap + conflict + tradelog 
                          + europe +africa + asia + america, model="negbin", 
           data=mi(imp1,imp2,imp3,imp4,imp5,imp6,imp7,imp8,imp9,imp10))
z.lwd.out <- zelig(terrorinclead~transferslog + govleft + democ + poplog + govcap + conflict + tradelog 
                          + europe +africa + asia + america, model="negbin", 
           data=reasonabledata)

x.low <- setx(z.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)-sd(reasonabledata$tradelog,na.rm=TRUE))
x.high <- setx(z.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)+sd(reasonabledata$tradelog,na.rm=TRUE))

x.lwd.low <- setx(z.lwd.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)-sd(reasonabledata$tradelog,na.rm=TRUE))
x.lwd.high <- setx(z.lwd.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)+sd(reasonabledata$tradelog,na.rm=TRUE))

s.out <- sim(z.out,x=x.low,x1=x.high)
s.lwd.out <- sim(z.lwd.out,x=x.lwd.low,x1=x.lwd.high)

junk<-sort(s.out$qi$fd)
ci5<-c(junk[round(length(junk)*0.05)],junk[round(length(junk)*0.95)])
m[5]<-mean(junk)

junk<-sort(s.lwd.out$qi$fd)
cilwd5<-c(junk[round(length(junk)*0.05)],junk[round(length(junk)*0.95)])
mlwd[5]<-mean(junk)

# Model 6

z.out <- zelig(terrorinclead~welfarelog + govleft + democ + poplog + govcap + conflict + tradelog 
                          + europe +africa + asia + america, model="negbin", 
           data=mi(imp1,imp2,imp3,imp4,imp5,imp6,imp7,imp8,imp9,imp10))
z.lwd.out <- zelig(terrorinclead~welfarelog + govleft + democ + poplog + govcap + conflict + tradelog 
                          + europe +africa + asia + america, model="negbin", 
           data=reasonabledata)

x.low <- setx(z.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)-sd(reasonabledata$tradelog,na.rm=TRUE))
x.high <- setx(z.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)+sd(reasonabledata$tradelog,na.rm=TRUE))

x.lwd.low <- setx(z.lwd.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)-sd(reasonabledata$tradelog,na.rm=TRUE))
x.lwd.high <- setx(z.lwd.out,tradelog=mean(reasonabledata$tradelog,na.rm=TRUE)+sd(reasonabledata$tradelog,na.rm=TRUE))

s.out <- sim(z.out,x=x.low,x1=x.high)
s.lwd.out <- sim(z.lwd.out,x=x.lwd.low,x1=x.lwd.high)

junk<-sort(s.out$qi$fd)
ci6<-c(junk[round(length(junk)*0.05)],junk[round(length(junk)*0.95)])
m[6]<-mean(junk)

junk<-sort(s.lwd.out$qi$fd)
cilwd6<-c(junk[round(length(junk)*0.05)],junk[round(length(junk)*0.95)])
mlwd[6]<-mean(junk)

ci.all<-cbind(ci1,ci2,ci3,ci4,ci5,ci6)
cilwd.all<-cbind(cilwd1,cilwd2,cilwd3,cilwd4,cilwd5,cilwd6)

my.main="First Differences of Trade Dependence on Violence"
my.xlim<-c(-1.5,1.5)
my.ylim<-c(-1.5,1.5)

plot(m,mlwd,xlim=my.xlim,ylim=my.ylim,main=my.main,xlab="imputed dataset",ylab="listwise deleted dataset",col="red")
#grid.rect(x = unit(1, "npc"), y = unit(-1, "npc"),  width = unit(1, "npc"), height = unit(1, "npc"), gp = gpar(fill = "light grey"),draw=TRUE)
rect(0,-3,3,0,col="light grey")
rect(-3,0,0,3,col="light grey")

par(new=TRUE)
plot(m,mlwd,xlim=my.xlim,ylim=my.ylim,main=my.main,xlab="imputed dataset",ylab="listwise deleted dataset",col="red")
par(new=TRUE)

for(i in 7){
for(j in 1:6){
  par(new=TRUE)
  plot(rbind(m[j],m[j]),cilwd.all[,j],xlim=my.xlim,ylim=my.ylim,main=my.main,xlab="imputed dataset",ylab="listwise deleted dataset",type="l",col="red")
  par(new=TRUE)
  plot(ci.all[,j],rbind(mlwd[j],mlwd[j]),xlim=my.xlim,ylim=my.ylim,main=my.main,xlab="imputed dataset",ylab="listwise deleted dataset",type="l",col="red")
}
}

abline(a=0,b=1)
abline(h=0)
abline(v=0)
 
